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Abstract 

Computational aeroacoustics requires efficient, high-resolution simulation tools. And for 
smooth problems, this is best accomplished with very high order in space and time methods 
on small stencils. But the complexity of highly accurate numerical methods can inhibit 
their practical application, especially in irregular geometries. This complexity is reduced by 
using a special form of Hermite divided-difference spatial interpolation on Cartesian grids, 
and a Cauchy-Kowalewski recursion procedure for time advancement. In addition, a stencil 
constraint tree reduces the complexity of interpolating grid points that are located near wall 
boundaries. These procedures are used to automatically develop and implement very high 
order methods (> 15) for solving the linearized Euler equations that can achieve less than one 
grid point per wavelength resolution away from boundaries by including spatial derivatives 
of the primitive variables at each grid point. The accuracy of stable surface treatments is 
currently limited to \\ th order for grid aligned boundaries and to 2 nd order for irregular 
boundaries. 


I. Introduction 

Noise generation and propagation are difficult to simulate numerically for a variety of 
well documented reasons, and require high-order numerical schemes. 1,2 However, high-order 
schemes can introduce a number of complications, such as: 

1. Large stencils near boundaries, with either many fictitious grid points, or large one- 
sided stencils, introduce programming complexity and numerical instability. 3 ’ 4 
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2. High order finite difference equations can require boundary treatments beyond the 
physical conditions of the original problem, which can excite spurious waves and in- 
stabilities. 0 

3. Generation of high-order, smooth, body-fitted grids around complex configurations can 
be difficult. 6 

4. High-order formulations can lack nonlinear robustness. 6 

5. The general usefulness of high-order methods is limited by first order accurate shock 
capturing. 7 

The approach described in this paper will focus on the first three issues. The fourth issue 
is being investigated. 8 And the last issue may possibly be avoided by including the physical 
viscosity and resolving the steep gradients, 9 but with very high resolution methods. 

In addressing the first three concerns, we limit ourselves to wave propagation and scat- 
tering problems governed by the linearized Euler equations. Previous studies point out the 
advantages of high order methods for acoustical propagation. 10 ” 15 Many methods in general 
use stop at 4 th order accuracy for time dependent problems since they use Runge-Kutta 
methods. High-order Runge-Kutta methods become notoriously difficult to derive because 
the number of nonlinear order conditions that need to be solved grows exponentially (i.e., a 
12 </l order method has 7813 nonlinear order conditions). 16 ” 21 The advantages of using Runge- 
Kutta methods at orders less than 6 are commonly cited as flexibility, large stability limits, 
and ease of programming. 22 The practical limit on their order has been an impediment to 
the analysis of their use in high order approaches for time dependent acoustic applications. 23 

In this paper we use a series of explicit , local, high order methods which have the same 
order of accuracy in space as in time. 24,25 These methods use Hermite interpolation on sten- 
cils that are two points wide, and a Cauchy-Ivowalewski recursive procedure 26 for obtaining 
time derivatives from the space derivatives of the interpolant. The time derivatives are then 
used to advance the primitive variables and their spatial derivatives in time with a Taylor 
series expansion. This general approach is called the Modified Expansion Solution Approx- 
imation (MESA) method. 11 This method can be used to derive and implement algorithms 
with arbitrarily high orders of accuracy in multiple space dimensions if their complexity is 
properly managed and the computer’s floating point precision is sufficiently high. 27 

The complex task of developing and coding a multidimensional interpolant for each MESA 
scheme can be eliminated using the tensor product 28 of a new divided-difference form of 
Hermitian spatial interpolation on a two point stencil. 27 And the task of obtaining time 
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derivatives with Cauchy-Kowalewski recursion 26 can be implemented by unrolling the recur- 
sion. 29 With these techniques, very high order MESA schemes may be implemented using 
only a few pages of code and their accuracy may be adjusted by merely changing loop indices, 
not the code. 27 This in turn enables local solution adaptive order changes at each time step. 

Most unsteady flow simulators using Cartesian grids in complex geometries are based 
on a finite volume approach 30 with the exception of the work by Kurbatskii and Tam, 31 
which specifically investigates acoustic scattering using a 4 th order accurate in time finite 
difference approach. In their work, curved surfaces are approximated with linear segments 
and ghost points are used to enforce boundary conditions; All of which incurred a fair amount 
of programming complexity to implement. 31 

The approach presented here, however, does not need to approximate the geometry and 
the programming task is completed automatically using computer algebra 29 in a manner sim- 
ilar to the works of Wirth 32 and Steinberg. 33 This is accomplished by viewing the boundary 
conditions as applying uniformly in time throughout the boundary surface. The boundary 
conditions can be differentiated in the surface and in time, and the governing equations can 
be used to obtain an infinite number of constraints in the boundary. 4 And a stencil con- 
straint tree is used to simplify the task of symbolically imposing these high order surface 
boundary conditions. 29 

With these techniques, the accuracy of interior propagation seems limited only by the 
floating point precision available. However, we stably attained only eleventh order accuracy 
for acoustic scattering in a box with sides that are aligned to the grid, and we attained only 
second order accuracy in more general cases. These limits are due to poorly conditioned 
matrix systems and numerical instability. However, these effects are dependent upon the 
choice of boundary conditions and require further study as will be shown. 

The objective of this paper is to present and validate a new approach to aeroaeoustic 
computations in complex geometries that has the potential to fully utilize the precision of 
today’s computers. A common theme in this paper will be the many advantages of using 
a two-point vide stencil which include very high resolution, solution adaptability, and CFL 
stability near boundaries. One of the more intriguing results to be presented here is the 
possibility of subgrid scale resolution using solution adaptable algorithms. 

II. Governing Equations 

We will demonstrate this new approach by solving the linearized Euler equations in a 
uniform mean flow field with mean velocity, M, perturbation velocity, u, pressure, p, and 
Cartesian coordinates x,. We assume the initial conditions are known and that no additional 
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sources are present. The conservation of momentum and energy equations are: 


dt 1 dx i dxj 

dp A J dp duj 

dt. ‘dx i dxj 


= 0 
= 0 


( 1 ) 


and are nondimensionalized with respect to the following scales, length scale L, velocity scale 
c (sound speed), time scale L/c, density scale p 0 ( ambient density), pressure scale p 0 c 2 . Also 
note that the continuity equation is not included since it is not necessary for calculating the 
acoustic response. 

We ’null use 2s + 1 order explicit MESA methods, where s can be any nonnegative integer. 
We will use only two-point stencils, where each grid point contains the primitive solution 
variables p and u, and their spatial derivatives (Hermitian data), for a total of 3(s + l) 2 or 
4(.s- + l) 3 terms per grid point in two or three spatial dimensions, respectively. 


III. Grid Classification 

Generally, we are interested in simulating the acoustical scattering from objects that arc 
defined parametrically, and which are superimposed on a uniform Cartesian grid, as shown 
in Fig. 1 for a cascade of three airfoils. Notice each grid point is labeled with a closed circle, 
an open circle, or the letter ”B” and represent an interior, a. fill, or a boundary grid point 
respectively. In the figure the assumed stencil size is three points per dimension and a fill 
grid point is simply one in which one of its neighboring grid points is either on or beyond 
the boundary. It is referred to as a fill since that grid point cannot be computed directly and 
need* to be '■filled” with data. If none of the neighboring grid points is on or beyond the 
boundary then it is classified as an interior grid point and can be computed normally. All 
other grid points are considered boundary points and are not needed since ghost grid points 
an* not used. 

Grid point classification proceeds by first labeling all grid points as boundary type, and 
then finding the interior and fill type grid points with a simple recursive procedure. The 
procedure determines if any neighboring grid points are boundary points by checking if the 
surface of any object intersects the imaginary line joining the center point to its neigh- 
bor. Since we are using a finite-difference method we need to check only for line-to-surface 
intersections, which is much simpler than determining the surface-to-surface intersections 
necessary if a finite-volume method is used. If there is no intersection, then the current grid 
point is a fill point and recursion stops. Otherwise it is an interior point and the procedure 
is called again, but starting with the neighboring grid points. 
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IV. Advancing the Solution using MESA 

With the uniform Cartesian grid points correctly labeled and the initial data assigned to 
all the interior and fill grid points, the next step is to advance the solution in time. Since 
interior grid points and their stencils are completely contained within the computational do- 
main, they will be advanced directly using an efficient form of the MESA schemes. However, 
since a fill point’s stencil reaches outside the computational domain field, its data will be 
obtained by spatial interpolation. At each time step, every grid point is first treated as an 
interior point to simplify the looping structure for efficient vector and parallel operations, 
and then the relatively few fill points are recomputed using interpolation. 

Advancing the interior points with MESA is a three step process. The spatial deriva- 
tives at the center of the stencil are approximated using spatial interpolation, the spatial 
derivatives are converted to time derivatives, and the solution is advanced in time with a 
Taylor series. The MESA scheme 34 is a general approach to developing high order numerical 
schemes, but for this work we are only interested in methods with stencils that are two 
grid points wide in all space dimensions because they enable arbitrarily high order spatial 
interpolation with no decrease in CFL stability bounds near a boundary. 

These two-point stencil schemes are numerically stable only when a staggered grid is 
used at each time step 34 as depicted in Fig. 2 for a two-dimensional stencil. The locations 
marked with ”X”s in the figure are advanced first using the 2x2 stencils centered about 
each ”X”. Next, the Hermitian data now stored at the ”X” locations are used as initial data 
to advance the solution to the center as indicated by the large dot in the figure. Staggering 
the grid has the same effect as applying the MESA scheme to the entire 3x3 stencil in the 
figure and then adding artificial dissipation. It is however, more efficient to stagger the grid 
because neighboring grid points will reuse the data at the ”X” locations in the figure. In 
two-dimensions this staggered grid procedure is numerically stable if: 34 


At 1 

Ax - 1 + max{|M x |, |A/j,|}‘ 


(2) 


A. Hermitian Derivative Approximation 

The MESA scheme requires an approximation to the solution of the primitive variables 
(p and Ui) and their spatial derivatives at the center of each stencil. Let the function, f(x,y), 
be an approximation of one of the primitive variables from Eqs. (1) and define its origin 
to be at the center of each two-point stencil shown in Fig. 3. Then in two-dimensions the 
following data must be approximated once for each primitive variable (f(x,y)=p(x,y), u(x,y), 
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and v(x,y)): 


d" +fe /( 0,0) 

dx a dy b 


V a, b : a, 6 


0, 1, 2, . . . ,2s + 1 


(3) 


On a Cartesian grid these approximations can be found using tensor product interpola- 
tion 28 in two or three dimensions if the interpolant: 


2s+l 2 s+1 

ffav) = £ H (-1) 

2—0 j = 0 

with its origin at stencil center is used with a 2s -1- 1 order method and if the following data 
is available at each grid point: 


d"+ fe /(.-r, y) 
dx a dy b 


V a, b : a, 6 = 0, 1 , 2 , . . . , s 


( 5 ) 


Determining the C{ 3 = terms in the two-dimensional interpolant (Eq. 4) 

for very high order methods is very difficult and inefficient unless it is reduced to a series 
of one-dimensional interpolations by using tensor products. Briefly, this is accomplished by 
performing one-dimensional interpolations of the form: 2 ' 

2s+l 

/M = E C’x' (6) 

2=0 


It is then possible to interpolate in the x direction, using the data on the grid points, the 
following two sets of data for y =y 0 and y = yp 


d i+j f{x = 0, y) % = 0, 1, 2, . . . , 2s + 1 
dx'dyi ' 3 — • ■ ■ ) s 


(7) 


Then this data is used to interpolate in the y-direction using the following one-dimensional 
shape function: 

2s+l 

f{x = 0,y) = Y, C jy J ( 8 ) 


Once the Cj terms arc determined, they provide the data shown in Eq. (3): 



/ 1 \ d'+i f(x = 0, y = 0) 

yd.j! J dx'dyi 


(9) 
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The one-dimensional shape functions (Eq. (6) and Eq. (8)) can be directly solved using 
computer algebra (Groebner Bases) 18 to create algebraic expressions for each term for up 
to approximately 30 ,; ' order methods. However, this results in lengthy, inefficient expressions 
that limit the accuracy and prevent instantaneous accuracy changes necessary for resolving 
varying wavelength scales. A better way is to use Hermitian divided-differences which will 
create the equivalent one-dimensional shape function (Newton’s interpolator}’ formula): 35 

f( x ) = QiA x ~ x oY +1 { x ~ x iY~ (s+l) 

+ '£,QiA x - x oY 0 °) 

i = 0 

where the are the traditional Hermitian divided differences from the tableau (shown in 
Fig. 4) which are given by: 


Do[Do[ 


n - (±\ #/ (g o) . 

1 0 ~ \j< ) dxi ’ 


Qi+s+ IJ 

> {h o, *}] 

, {/, 0, s}], 


(jr) 


dxl ’ 


Do[Do[ 

(~\ (Q tj — i ~ Qj - 1 j - * ) 

— Ax 

, {j, i - s, •;}] 

, {i, s + 1,2 * s + 1}]; 


Derivatives Algebraically Evaluated 


It is not enough to simply find Newton’s interpolator}’ formula, Eq. (10), since we need 
to approximate p, u, and v and their spatial derivatives at the stencil center. This would 
ordinarily require differentiating Newton’s Eq. (10) using computer algebra. By the product 
rule of differentiation, Newton’s formula will double in size after each differentiation, and 
a 50 th order MESA scheme would need a Newton’s formula of approximately 2 50 ~ 10 16 
terms. This issue is eliminated by the following result which is applicable only to two-point. 
stencils: 27 


d*7(g = 0) 

dx dx 


2s+l 

^ . Qi,i Ei y d x, 

i—dx 


(ii) 
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where the function Z itdx is independent of space and time, and can be computed simply as 
follows: 


Do[Do[ 

Zi,dx — 


i\ 


,{t,<ic,s}],{<fa.0,(2*+l)}] 

Do[Do[ 

%i4* — 

dx J r ] 

g 0 ( ( tfa-r)!r ! (- 3:(, ) (,41 - r> (- Il ) (i -‘- 1 "' ,+ * 

n£; , -'[t-(« + i)-e]* 

IKo [»+!-*]) 

, {t, s + 1,2s + 1}], {dx, 0, (2s + 1)}] 


( 12 ) 


(13) 


With these developments it is now possible to efficiently approximate all solution variables 
and their derivatives at the center of a two-point stencil to any level of accuracy with only a 
page of code. Since these explicit forms are completely algebraic it is possible to dynamically 
adapt the accuracy of the approximation to the solution evolution. 

B. Time Evolution 

To fully utilize the arbitrarily high order accuracy in space that is possible on two-point 
stencils, it is necessary to achieve similar accuracy in time; This can be accomplished with the 
MESA method. This method uses the governing equations, Eqs. (1), to convert the recently 
approximated spatial derivatives to mixed space-time derivatives in a manner reminiscent of 
the Lax-Wendroff approach 36 by differentiating the following equations in space and time: 



M f du i , 

= - (AI - aJT + 


dp 

Ot 


< M ‘^ + 


dp 

dxj 

duj 

dxj 


(14) 


The mixed space-time derivatives are then used in modified series expansions with local 
coordinates about the center of the interpolation stencil at the current time level. For two- 
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dimensional problems, the following series are used, with O = 2s + 1: 


V, t) 
u(x,y,t) 
v{x, y, t) 


O O 2 ( 0 ) 

EIE c?>y 

2—0 j — 0 k ~ 0 
O O 2 ( 0 ) 

£ £ £ cy^Vt* 

2=0 j=0 A~=0 
O O 2(0) 

£ £ £ 

2=0 J=0 k~0 


( 15 ) 


„.l irrf . rf = i 1 1 &*' t 9 + */(a=0.,=0,t=0) 

^ z j , A V i !j ! Ar ! ' &T 1 Of* 1 ’ 

The variables p, u, and v, and their spatial derivatives, are then advanced to the next 
time step by evaluating Eq. (15) and its spatial derivatives with (x = 0, y = 0, t — At), as 
follows: 


d fl+6 p(0,0, At) 

dx a y b 

d a+6 u(0, 0, At) 

dx a y b 

<T ,+fe t>(0,0, AQ 

dx a y b 


2 ( 0 ) 

= £ (oWJCJ^Af* 


jfc=0 

2 ( 0 ) 


= £ (a'W)C« Ai ‘ 


0 

2 ( 0 ) 


= £ (aWX^Al* 


(16) 


By unrolling the Cauchy-Kowalewski recursion in Eq. 14, we can quickly express the 
mixed space-time derivatives C% bk , C^ bk , and C^ b k in Eq. (16) in terms of the space deriva- 
tives C % b 0 , Q m , and Q 6>0 , which were approximated at the stencil center in the last section. 
The following loop efficiently unrolls the recursion: 29 


Do[Do[Do[ 

a = a + l]b = b + l] k = k — 1 

(~a,b,k = 

(-((a),(M,.Cf kt + C» a )) 


(*) * (My * c :.U + C li^' k 


(17) 
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c u — 

L a,h,k ~ 


(-((V) * My * C™ j j) - 

> 

'a,M J 


(«) * (q a + w. ♦ cj Ai ))/* 


y^L" 

L 'a, 6,* ~~ 

(-«»)*(<=:«+ a/, *cj w )) 
(«) » m, » q A( )/fc 

o,0,0],6,0,0],t,l,2»0]; 


This simple form applies to uniform mean flow, but for general flows Cauchy-Kowalewski 
recursion becomes complicated. 26 For example, using Eq. (14) to calculate the mixed space- 
time derivative, C“ 01 , when the mean flow varies in space, results in the expression: 


d 2 u dMi du , , d 2 u <9 2 p. 

dtdx dx dxi + 1 1 dxf + dx 2 


(18) 


which requires knowing the derivatives of the mean flow, and higher order derivatives of 
Eq. (18) require higher order derivatives of the mean flow as well. Getting this information 
may require computing the mean flow using a very high order MESA method as well and 
the equations will grow exponentially unless a simple form analogous to Eq. (11) is found. 
Despite these complications, Cauchy-Kowalewski recursion, when automated with computer 
algebra tools, 29 can produce numerical schemes with higher order than is currently possible 
with Runge-Kutta methods. 


V. Including Wall Boundaries 

The procedure for advancing the interior grid points is not applicable to fill grid points 
near boundaries since part of their stencil is on or within a solid object. We obtain data 
for these points at each time step with a Hermite interpolant that is consistent with the 
wall boundary conditions, and which uses the nearby solution at interior grid points. This 
Hermite interpolant is equivalent to the interpolant used for the interior grid points in Eq. (4), 
but is written as: 


f{x,y) = 


1 5 


£ £ Wv)) 

dx^dy — 0 


d ix+dv f{xi,y j ) 
dx dx dydy 5 


(19) 


where the H Xudx {x) and H yj j y (y) terms are 2s + 1 order polynomials in x and y respectively. 
Since the data at the interior grid points is known, this Hermite form of the interpolant has 


NAS A/TM— 2000-2 10378 


10 


the advantage of reducing the number of unknown coefficients, Cfj, in Eq. (4), to only the 
unknown data which is required at the fill points, - 

Each shaded box in Fig. 5 represents an interpolation region used for interpolating the 
values at the fill points in each box. The boundary conditions are imposed upon the shape 
function for each box at the locations on the surface intersected by the arrows. The numbers 
inside each box show the sequence in which the fill points are interpolated. 

A. Choice of Wall Boundary Conditions 

Each fill point on the grid provides (s + l) 2 unknowns in Eq. (19), so that an equivalent 
number of constraint equations must be obtained. An infinite number of constraints can be 
obtained by exploiting the fact that the physical boundary conditions apply uniformly in 
time everywhere in the boundary. 4 For Eqs. (1), the inviscid boundary condition is that the 
normal component of velocity is zero at a surface, 

V-fj = 0. (20) 


Since vorticity is convected with the mean flow, we also assume in this work that there is no 
vorticity at a surface, 


dv du 
U dx dy 


For flat walls, the boundary conditions for this case are: 


( 21 ) 


Q2n+\ +T p 

dfj 7n+l df T 

d 2n+T Yr, 

df) 2 n df T 

d 2n+l+Ty f 

df] 2n+l df T 


= 0 

= 0 , 
= 0 


Vrc, T : 

n,T e (0,1,2,...) 


( 22 ) 


where V is the magnitude of the perturbation velocity vector with velocity components u and 
v in the Cartesian x arid y-axis directions, p is the scalar pressure, fj is the unit vector normal 
to the wall surface direction, f is the unit vector tangent to the wall boundary direction, If 
is the velocity tangential to the wall, and is the velocity normal to the wall. 

Selecting which set of boundary conditions to use and where to apply them on the 
boundary can be simplified by considering the general form of the normal and tangential 
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(23) 


derivatives of f(x,y): 37 

d N+T f(x, y) _ 
df] N dr T 

(Vx£ + %%) N {-Vy£ + Vx^) T f{x, y ) 

If each operator is expressed in its binomial series form and the order of summation rear- 
ranged, then by substituting Eq. (19) into Eq. (23) the following form is derived: 


d v+7 7 (x,y) 

QfjNfT 


£ £ 

0 dx ? dy — 0 


[n 




&**+*« f(x u ?y .) 




dx,dy^Xi ,yj J Qx^ x dl)^ 


(24) 


where Y dx , d y.x„ yj is defined by: 


\ r 

1 dx,dy,Xi,yj 

££*■(«, b) 

a— 0 6—0 

K(a,b) = 

N !T! 


d N * T (H,„ dI (x)B a ,, dy (y)) 


Q ^ A r -f-T— a—b 6 


(25) 


(26) 


f A 7 IT! \ 

\ (A 7 — 6)!M(T— a)!a! ) '>■ 


N + a - b , n T— a+b 


'y 


'(-i) 


T-a 


Since maintaining numerical stability for hyperbolic problems normally requires methods 
that include all information from within the domain of dependence defined by the character- 

QN +T {jj ^ d ) 

istic surfaces, 38 we want to select N = 2n + l and T so that 0x N + rYa-h d yi+i in Eq. (25) never 
becomes zero, because this would exclude some of the influence of neighboring grid points. 
This is accomplished by choosing (n : n = 0,1,2, ..., s) and (T : T = 0, 1, 2, 2(s — n)) for 
the first and third boundary conditions in Eq. (22). Those conditions will produce linearly 
consistent systems that can be solved to determine the interpolant, or equivalently, the data 
required at the fill point. However, in the special case in which rj x = 0 or r) y = 0, we need 
to use the different conditions (n : n = 0, 1 , 2, ..., s) and (T : T = 0, 1, 2, ..., s) to insure non- 
singular matrix systems. This special case does not result in a loss of information because 
the affected terms in Eq. (25) are supposed to be zero. For example, a third order MESA 
method would use the following pressure boundary conditions at each location indicated by 
the arrows in Fig. 5: 


dp 

df] 


7TT = 0, 


d 2 p 


T = 0, 


dp 


dfjdf ’ dpdr 2 


-.g- 


(27) 
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The same method applied to the special case shown in Fig. 6 would use the conditions: 


dp _ 0 d 2 p _ 5 3 p _ d A p 

Of) 'dfjdf ’ drf ’ diff 

B. Symbolic Solution 

The shape function for each fill point can now be solved using the boundary conditions 
and the known interior data in each shaded region of Fig. 5. The boundary condition 
Eqs. (22) are applied to the shape function Eq. (24) to solve for the unknown data at the 
fill points. This linear system of (s + l) 2 equations per fill point may be written as: 

Mf = Nd (29) 

where M and N are matrices whose elements are given by Yd x ,dy,x t ,yj in Eq. (25), where / is 
a vector of the data needed for the fill points, and where d is a vector of known data from 
the interior points. 

This system can be solved using any linear solver at each time step, but for high order 
methods this becomes expensive and the matrices become ill-conditioned. A better approach 
is to symbolically solve this system once to form the algebraic solution: 

/ = M -1 Nd (30) 

where M and N are numerical matrices dependent upon the geometry but the vectors are 
symbolic. This results in an algebraic expression for each fill point as a linear combination of 
the interior grid point data. It is then a simple matter to evaluate the fill points by updating 
the vector d at each time step and evaluating their linear combinations. 

This approach worked well for up to lE h order accuracy in two dimensions for the case 
(shown in Fig. 6) of no mean convection in an unrotated box that is aligned with the co- 
ordinate axes. For higher order algorithms, the matrix M becomes ill-conditioned, and 
finding its numerical inverse becomes difficult, inefficient, and unstable. Resorting to Gaus- 
sian elimination would be more stable for high order cases, but the cost at each time step is 
prohibitive. 



VI. Stencil Constraint Tree 

Instead of using a poorly conditioned numeric matrix M, its symbolic form can be in- 
verted. There are actually only a relatively few symbolic matrices which need to be inverted 
regardless of the geometry! These few cases can be found by using a new tree data struc- 
ture, referred to as a Stencil Constraint Tree (SCT), which represents all possible stencil 
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configurations. For example, in Fig. 7, all possible 3x3 stencil configurations are shown in 
which: 

• the top, center grid point is a ’’boundary” point (indicated by ”B”), 

• the center, center grid point is a ’’fill” point (hollow circles), 

• and the bottom, center grid point is an ’’interior” point (filled in circles). 

The stencil constraint tree (SCT) for the simple case of a 2 x 2 stencil which has an interior 
grid point in its top-right location is shown in Fig. 8. It is constructed by propagating 
symbolic constraints in a manner analogous to Waltz’s symbolic constraint propagation 
algorithm. 39 

Once the tree is constructed, a particular stencil configuration is found by traversing the 
SCT from top to bottom along any path and converting the branch labels using Figs. 9 
and 10. For example, the branch number 10 corresponds to position 4 and label 1 in Fig. 9 
and represents an ’’interior” grid point that is located at the top right corner of the stencil. A 
similar procedure is applicable to larger stencils as well. Notice that identifying a particular 
n-point stencil using a SCT can be done in n 2 or n 3 steps for two and three dimensional 
stencils respectively instead of the 3 7,2 n 2 and 3" 3 n 3 comparisons typically required using a 
brute-force comparison of all stencil configurations. 

A. Constructing the Tree 

The most significant advantage of the SCT is it can reduce the number of symbolic 
matrices from Eq. (29) that need to be considered. This is accomplished by propagating 
natural constraints during the construction of the tree. 29 Natural constraints are simply a 
list of rules which are a natural consequence of the definitions of interior, fill, and boundary 
grid points on a Cartesian grid. For example, some natural constraints which limit the set 
of stencil configurations are: 

1. An interior grid point is never adjacent to a boundary grid point and its converse, 

2. And for two and three-point stencils, a fill grid point will always be adjacent to at least 
one boundary point. 

If the natural constraints are not enforced then the list of stencil configurations in Fig. 7 
grows to 3 6 = 729 since only the middle column is defined. The advantages of enforcing 
natural constraints increase with both spatial dimension and stencil size. 
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B. Small Stencils Guarantee CFL Stability 

These techniques can be applied to both of the 5x5 stencils shown in Fig. 11 to construct 
their SCT. This SCT will represent all possible 5x5 stencils which can occur with a fill point 
at its center. And we find that there arc only 65 unique rotationally symmetric cases of 3 x 3 
stencils centered about the fill point; 29 All other possible stencil configurations are simply 
rotations of this core group! Note that we assume there will be at least one interior point 
next to the fill point, otherwise that fill point is not needed as shown in the corners of Fig. 5. 
As shown in Fig. 11, there will always be a 3 x 3 stencil contained within the 5x5 stencil 
that does not contain a boundary point and yet always contains a fill point on its edge or 
corner. This guarantees the existence of a 3 x 3 stencil for each fill point in any geometry that 
is not intersected by the surface geometry! Therefore, it will always be possible to expand 
the stencil outward from each fill point to the surface where the boundary conditions are 
imposed. This in turn guarantees that the CFL stability criterion will be satisfied, since the 
domain of dependence for the fill interpolaiit is increased from a standard stencil size. 

This CFL guarantee is a property of two-point and three-point stencils in two or three 
dimensions. Four-point stencils or larger will require smaller time steps near the boundary 
and make it difficult to choose locations on the surface that guarantee a linearly consistent 
system in Eq. (29). 


VII. Mapping 

All possible 5x5 stencil configurations containing a fill point are known after the con- 
struction of the SCT for the stencils in Fig. 11. This information is used to select locations 
on rln- boundary at which the boundary conditions are enforced. A surface will always occur 
between a fill point and a boundary point, so that even before we know exactly where the 
surface is. we know in which direction to go to find the surface. 

We would like to use the points on the surface that are closest to the interpolation 
stencil while simultaneously insuring that never more than N locations are collinear with 
the coordinate axes for each Appoint interpolation stencil. These conditions maximize the 
accuracy of the interpolant and insure that the matrix M in Eq. (29) is nonsingular. One 
such mapping is shown in Fig. 12 for the upper triangular part of the outlined S8 3 x 3 stencil 
in Fig. 1 1 . The degenerate cases occur because it is not guaranteed that a nearby surface 
will be found using the mapping shown. These cases are usually not an issue in practice and 
can be avoided by choosing alternate interpolation regions as in Fig. 13 or by constructing 
the SCT for a larger stencil around the fill point in Fig. 11 and then choose a mapping with 
this additional information. 

Implementing higher accuracy boundary conditions with this approach requires increasing 
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the depth of Hermitian data at each grid point and increasing the number of boundary 
conditions used from Eq. (22), but the stencil size used for interpolation is never changed. 
Therefore, we never need to remap the fill points to the boundary and the CFL stability is 
not affected. 

A. Efficient Retrieval of Fill Point Solutions 

Once all the stencil configurations are known and the direction in which to map each fill 
point to the surface is found, then the symbolic form is known for every matrix M and N 
in Eq. (29). Each of these systems may then be symbolically solved and their subsequent 
solutions may be stored in the leaves of the SCT for fast retrieval later. 

The leaves of the tree will contain algebraic solutions for every fill point in any geometrical 
configuration. They can be used and reused many times without the need for solving Eq. (29) 
again. For example, if the walls are moving then the Cartesian grid points will need to be 
relabeled, but the same SCT, mapping, and symbolic solutions can be used at each time 
step to efficiently interpolate the data needed at each fill point. 

VIII. Numerical Results 

All of the following results solve the linearized Euler equations, Eqs. (1), in uniform mean 
flow using the approaches previously described. This approach to developing very high order 
methods was applied to both wave propagation and acoustic scattering problems in various 
test geometries. 

A. Acoustic Propagation 

For mean convection velocity M in an open bi-periodic (tri-periodic) domain (for d=2 or 
3 dimensions), the linearized Euler equations have the following analytical solution: 

d d 

p(x, t) = C t n S Xi , Ui (x, t) = - StjC x , n S Xi (31) 

l— 1 j & 

J =1 


where 


C t = 
ITT 1 1 = 


cos(IFf7rfj), 

cos(7rf||TT||), 




S Xi = sin (TT7r.fi) 

„ IT sin(7r/||ir||) 

^ f 7 — t 

Ill'll 

ft = (xi - Mil) 


The properties of very high order methods were tested on this problem with mean ve- 
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locity M = (1,1) and wavenumber IT — (1,1), and the results were used to validate the 
automatically generated codes. The methods were compared for efficiency and resolution; 
the results for up to 29 th order accuracy (using 64-bit precision) are shown in Fig. 14. Notice 
that higher-order methods perform better in each figure until machine precision is reached. 
By increasing the wavenumber in Eq. (31) to IT = (TFi > 2, TF 2 > 2), we also see that very 
high order methods offer subgrid scale resolution as shown to the left of the ordinate axis at 
the top of Fig. 14. Similar results were found while solving three-dimensional cases as well. 29 

The resolution of the methods improves while using 128-bit precision as the accuracy 
increases until about the 2l s/ order for signals with 8 grid points per wavelength; This limit 
varies with signal frequency. If the order is increased above this limit, then the resolution 
of each signal degrades from roundoff errors v r hich begin to dominate around 57 th order 
accuracy for signals with a wavelength of 8 grid points. However, by selecting wavenumber 
IT — (16, 16) in Eq. (31), and by using 40 th — 6Q ih order methods on a two-point stencil, 
it is possible to resolve 4 wavelengths per grid point. Note that this is not a violation 
of the Nvquist criterion 41 since each grid point contains multiple data values for Hermite 
interpolation. 

For extremely high order methods, the convergence slope at the top of Fig. 14 is so 
steep that essentially only a specific frequency is correctly resolved wfith finite precision 
floating point hardware unless a method for detecting roundoff error is available. Fortunately, 
roundoff error may be detected by checking the columns of the Hermite divided-difference 
tableau in Fig. 4, w r hich is used for the interpolation. The sum of column j should equal the 
difference of the first and last terms in column j-1. 27,40 If roundoff error begins to occur, 
simply stop constructing the tableau at that point and use the resulting lower order method. 
This saves computational effort and enables the solution of a field containing widely varying 
wavelengths. 

B. Acoustic Scattering 

With the exceptional qualities of very high order methods established for wave prop- 
agation, the simulation of acoustic scattering from within a rotated box and a circle was 
completed as showm in Fig. 15. An analytical solution to wave scattering within the rotated 
box with no mean convection is: 


p{x, V , t ) 

u(x, y, t ) 

v{x,y,t) 


— COs(v^7T t) COs(n x) cos(7r y) 
cos(7r y) sin(>/2 7 rt) sin(7rx) 

V2 

cos(7rx) sin(v^Trf) sin(7ry) 

72 


(32) 

(33) 

(34) 
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where the coordinate system is aligned with the box. An analytical solution to wave scat- 
tering within a circle with no mean convection is: 


p{r, 0 , t) 


5J 0 (Ar) cos(Af) 
v/’T-J o(A) 


(35) 


where A = 3.83171, J 0 is the Bessel function of the first kind of order 0, and where polar 
coordinates are used. 

For the case of an unrotated box we found stable solutions with second, third, fifth, 
seventh, ninth and eleventh order accuracy using boundary conditions from Eq. (22) with 
(n ; n = 0,1,2,...,$) and (T : T = 0,1,2,...,$) as shown in Fig. 17. Using the boundary 
conditions in Eq. (22) with n,T — 0 we found stable 2 nd order solutions for all box and 
circle cases, regardless of the rotation angle. The first five plots in Fig. 16 correspond to box 
rotation angles: 0, respectively. The last curve in the figure shows a 2 nd order 

method without wall boundaries. 


C. Numerical Stability 

In Fig. 18, both rotated boxes are identically solved with a 2 nd order method. The 
sequence in which the fill points are interpolated in each box is switched as indicated by the 
numbers in the center of each shaded region. The sequence used on the top box in the figure 
is numerically stable while the sequence used for the box on the bottom is not. It was found 
that hir 2 ,ld order methods, ordering regions by the number of fill points they contain from 
fewest to most and then interpolating the fill points in that order was always stable for the 
cases Tested. This is why most of the interior is shaded for the box cases in Fig. 15 compared 
to the boxes in Fig. 18. But, clearly, there are multiple stable choices of shaded regions and 
interpolation sequences. 

Rather than relying on numerical experiments, an important next step is to use the 
symbolic form in Eq. (29), of which only a relatively few cases are possible, to perform a 
complete stability analysis so that higher order conditions can be made reliably stable. 


IX. Conclusion 

Small stencils containing Hermitian data were found to possess key properties necessary 
for developing exceptionally efficient, explicit algorithms of arbitrary formal accuracy in space 
and time for scattering problems with arbitrary surface orientations. This is accomplished 
by approximating spatial derivatives using a special form of two-point Hermitian divided- 
difference spatial interpolation and an unrolled Cauchy-Kowalewski recursion procedure for 
time advancement. A stencil constraint tree is used to find all possible stencil configurations 
and their respectively well-posed boundary conditions so that grid points near arbitrarily 
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oriented surfaces may be interpolated. Computer algebra is then used to find the symbolic 
shape function for each stencil configuration, and the stencil constraint tree is used to quickly 
identify the correct interpolant for each near boundary grid point. 

In spite of the development of simple procedures for the time evolution at all grid points, 
the Fortran application code of a very high order method is complicated, so the task of 
programming was automated using computer algebra procedures. These procedures have 
successfully produced: (1) stable, parallel 57 th order wave propagation methods in two and 
three spatial dimensions with periodic boundary conditions; (2) stable 2 nd order methods for 
scattering problems with generalized surface boundary conditions in two spatial dimensions; 
and (3) stable \ V h order methods for scattering problems with surfaces aligned to the grid 
in two spatial dimensions. 

The resolution and efficiency of these methods improve with order until the accuracy 
exceeds the limits of machine precision. At these very high orders of design accuracy it is 
necessary to control roundoff error and a method for doing this is suggested. 

The procedures discussed here provide a systematic method for quickly changing the 
accuracy of both the interior propagation and wall boundary condition implementations to 
suit flow conditions. And it was shown that very high order methods (> 15) provide an 
opportunity for subgrid scale resolution. 

These results demonstrate the potential value of very high accuracy in time as well as 
in space for aeroacoustic calculations. But, a detailed stability analysis for these procedures 
remains to be done for the generalized high order surface conditions. And, a method for 
unrolling nonlinear Cauchy-Kowalewski recursion needs developed to achieve arbitrarily high 
accuracy in time for nonlinear problems. 

Nonetheless, it is necessary to use Hermitian data for achieving very high order solutions 
because Lagrangian stencils simply become too large. And since small stencils allow for a 
very simple interpolation procedure of arbitrary accuracy while maintaining CFL stability 
near boundaries, future work will be devoted to exploiting them in substantial applications. 
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Figure 6: Stable Unrotated Box Case 



NASA/TM— 2000-2 1 0378 


27 



Figure 7: All Possible S7 Stencil Configurations with Fill at Center : 16 Cases 



Figure 8: Two-Point Stencil Constraint Tree 
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Figure 10: Stencil Grid Positions 
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• = Interior Grid Point Type 
o= Fill Grid Point Type 
Z= Not a Boundary Grid Point Type 
X= Unspecified Grid Point Type 
B= Boundary Grid Point Type 



Figure 11: Sub 3x3 stencil symmetry 



Figure 12: S8 Symmetrical Mapping: 18 Cases 
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Figure 16: Resolution of Rotated Boxes and Circles 




Figure 17: Unrotated Box Grid Resolution Studies 
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